#include<stdio.h>
#include<stdlib.h>
#include<math.h>

extern __shared__ float x[];

/* Função que calcula k iterações do algoritmo na GPU */
__global__ void calcula(int tam, int num_ite, float *MA2, float *VB2, float *xfinal) {
        int i, j, k, num, tid = threadIdx.x;
        float soma = 0;

	/* primeiro setamos a primeira linha da "matriz" x, que aqui é um vetor */
        if(tid == 0){ 
                for(j = 0; j < tam; j++)
                        x[j] = VB2[j];
        }   

        __syncthreads();

	/* agora cada thread será responsável pelo calculo de novas iterações */
        if(tid<tam){
		num = (tam/500);
                for(k = 0; k < num_ite; k++){
			for(i = 0; i < num; i++){
		                soma = 0;
		                for(j = 0; j < tam; j++)
		                        soma += MA2[((((tid*num)+i)*tam)+j)] * x[j];
		                xfinal[((tid*num)+i)] = VB2[((tid*num)+i)] - soma;
			}

			if(k != (num_ite-1)){
				__syncthreads();

				for(i = 0; i < num; i++)
                        		x[((tid*num)+i)] = xfinal[((tid*num)+i)];

                        	__syncthreads();
			}

                }   
        }   
}

int main(){
        int i, j, k, J_ORDER, J_ROW_TEST, J_ITE_MAX;
        float J_ERROR, diag;
        float resA = 0, resB;
        float **MA, *VB, *MA2, *VB2;
        float *xfinal, *dev_MA2, *dev_VB2, *dev_xfinal;

        /* Leitura dos dados */
        scanf("%d", &J_ORDER);
        scanf("%d", &J_ROW_TEST);
        scanf("%f", &J_ERROR);
        scanf("%d", &J_ITE_MAX);

        /* Separa memÃ³ria dos vetores que utilizaremos, MA2 e xk sÃ£o enviados, xk1 Ã© retornado */
        cudaMalloc((void**)&dev_MA2, J_ORDER * J_ORDER * sizeof(float));
        cudaMalloc((void**)&dev_VB2, J_ORDER * sizeof(float));
        cudaMalloc((void**)&dev_xfinal, J_ORDER * sizeof(float));

        /* AlocaÃ§Ã£o de memÃ³ria das variaveis */
        MA = (float**)malloc(J_ORDER * sizeof(float*));
        MA2 = (float*)malloc(J_ORDER * J_ORDER * sizeof(float));
        for(i = 0; i < J_ORDER; i++)
                MA[i] = (float*)malloc(J_ORDER * sizeof(float));
	
        VB = (float*)malloc(J_ORDER * sizeof(float));
        VB2 = (float*)malloc(J_ORDER * sizeof(float));
        xfinal = (float*)malloc(J_ORDER * sizeof(float));

        /* Leitura da matriz MA */
        for(i = 0; i < J_ORDER; i++)
                for(j = 0; j < J_ORDER; j++)
                        scanf("%f", &MA[i][j]);

        /* Leitura do vetor VB */
        for(i = 0; i < J_ORDER; i++)
                scanf("%f", &VB[i]);

        /* Criacao de  VB* (neste caso VB2) */
        for(i = 0; i < J_ORDER; i++){
                diag = MA[i][i];
                VB2[i] = VB[i]/diag;
        }

        /* Criacao de MA* (neste caso MA2) */
        for(i = 0; i < J_ORDER; i++){
                diag = MA[i][i];
                k = 0;
                for(j = (i*J_ORDER); j < (i+1)*J_ORDER; j++)
                        MA2[j] = MA[i][k++]/diag;
                MA2[i*(J_ORDER+1)] = 0;
        }

        /* Copia todo MA2 e VB2 que serÃ£o utilizado */
        cudaMemcpy(dev_MA2, MA2, J_ORDER * J_ORDER * sizeof(float),cudaMemcpyHostToDevice);
        cudaMemcpy(dev_VB2, VB2, J_ORDER * sizeof(float),cudaMemcpyHostToDevice);

        calcula<<<1,500,(J_ORDER*sizeof(float))>>>(J_ORDER,J_ITE_MAX,dev_MA2,dev_VB2,dev_xfinal);

	/* Retorna xfinal */
        cudaMemcpy(xfinal, dev_xfinal, J_ORDER * sizeof(float),cudaMemcpyDeviceToHost);

        /* Calculo da resposta */
        for(i = 0; i < J_ORDER; i++)
                resA += MA[J_ROW_TEST][i] * xfinal[i];
        resB = VB[J_ROW_TEST];

        printf("Iterations: %d\n", J_ITE_MAX);
        printf("RowTest: %d => [%f] =? %f\n", J_ROW_TEST, resA, resB);

        /* Libera memÃ³ria CUDA */
        cudaFree(dev_MA2);
        cudaFree(dev_VB2);
        cudaFree(dev_xfinal);

        /* Libera memÃ³ria */
        free(VB);
        free(VB2);

        for (i = 0; i < J_ORDER; i++)
                free(MA[i]);
        free(MA);
        free(MA2);
        free(xfinal);

        return 0;
}

